Chaotic behavior in classical Yang-Mills dynamics 
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Understanding the underlying mechanisms causing rapid thermalization deduced for high-energy 
heavy ion collissons is still a challenge. To estimate the thermalization time, entropy growth for 
classical Yang-Mills theories is studied, based on the determination of Lyapunov exponents. Distinct 
regimes for short, medium and long sampling times are characterized by different properties of their 
spectrum of Lyapunov exponents. Clarifying the existence of these regimes and their implications for 
gauge-field dynamics is one of the results of this contribution. As a phenomenological application we 
conclude that for pure gauge theories with random initial conditions thermalization occures within 
few fm/c, an estimate which can be reduced by the inclusion of fermions, specific initial conditions 
etc. 



I. INTRODUCTION 

Experiments have shown that a new form of strongly 
interacting matter with very high energy density and un- 
usual transport properties is created in collisions between 
heavy nuclei at energies attainable at the Relativistic 
Heavy Ion Collider (RHICL up to 200 GeV per nucleoli 
pair in the center of mass [l[ . Theoretical arguments as 
well as circumstantial experimental evidence suggest that 
this matter is a strongly coupled quark-gluon plasma Q ■ 
The early thermalization of this matter leading to the for- 
mation of a quark-gluon plasma is one of the largest un- 
explained puzzles in RHIC physics. Hydrodynamic sim- 
ulations are consistent with a thermalization time of 1.5 
fm/c or less [||. It is generally believed that the instabil- 
ity and consequent exponential growth of intense gluon 
fields would be the origin of early thermalization. Vari- 
ous plasma instabilities such as the Weibel instability Q 
and the Nielsen-Olesen instability @ can cause the ex- 
ponential growth of the amplitude of unstable modes of 
the SU(3) gauge field. The plasma instability may be 
characterized by the negative curvature of the potential, 
leading to the equation of motion 



(1) 



where Xi denotes the field variable in the unstable mode. 
The energy stored in the intense gauge field eventually 
produces abundant particles and evolves towards a ther- 
malized state. The thermalization mechanism governing 
this transition is not yet clear, and the time scale on 
which it occurs is not known. The equilibration problem 
is simplified, however, by the high occupation probabil- 
ity of the unstable modes, which makes a quasi-classical 
treatment of the thermalization process, at least of its 
initial stages, possible. 

In the classical dynamics, the apparent entropy of an 
isolated system is produced by the increasing complexity 
in phase space. The distance between classical trajec- 
tories starting from very similar initial conditions grows 



exponentially in the long-time evolution of a chaotic sys- 
tem, 



\6Xi(t)\ oce 



(2) 



where 8Xi represents the separation of trajectories, and 
Xi is referred to as the Lyapunov exponent (LE). The en- 
tropy production rate is given by the Kolmogorov-Sinai 
(KS) entropy, which is defined as the sum of positive 
LEs, dS/dt = Sks = Y^x>o^i- The production of en- 
tropy at the quantum level poses additional problems 
such as the decoherence of the quantum state of the 
system Q, since the evolution in pure state generates 
no entropy and some kind of coarse graining is neces- 
sary. Kunihiro, Miiller, Ohnishi, and Schafer Q (hence- 
forth referred to as KMOS) proposed to apply the Husimi 
function, a smeared Wigner function with minimal wave 
packets, to define a minimally coarse grained entropy, the 
Wehrl entropy, and showed that it grows at the rate of 
the KS entropy in the classical long-time limit, i.e. if the 
system has enough time to sample the complete phase 
space. In that paper application was limited to simple 
cases where the number of degrees of freedom is essen- 
tially one. In this contribution we extend the KMOS 
framework to more realistic processes. 

In this work, we analyze the chaotic behavior in the 
classical Yang-Mills (CYM) evolution. Specifically, we 
analyze the exponentially growing behavior of the dis- 
tances between the trajectories. We find that we have 
to distinguish different regimes, depending on sampling 
time, namely a kinetic stage for short sampling times, an 
intermediate- and a long-time regiem.In each case we con- 
sider the exponential growth rate of the distance between 
two trajectories, which follows the equation of motion, 



SX(t) = H(t,X)6X(t) 



(3) 



where % is the so-called Hesse matrix or Hessian and 
analyse the time evolution of the distance vector 8X. 
i) The instantaneous change of 5X is determined by the 
eigenvalues of the Hessian, which we will refer to as the 



2 



local Lyapunov exponents (LLE). 

ii) The evolution of the distance on ergodic time scales is 
described by the Lyapunov exponents ([2]), which we will 
refer to as global Lyapunov exponents (GLE). 

iii) For the third, intermediate time period, the Hessian 
changes due to the nonlinear coupling among the different 
field modes, but the energy remains localized among the 
primary unstable modes. By using the Trotter formula, 
we can numerically integrate the equation of motion for 
the tangent space SX, and construct the time-evolution 
matrix for an intermediate time period. We will refer to 
the eigenvalues of the time-evolution matrix as interme- 
diate Lyapunov exponents (ILE). 

Because the ILEs describe the evolution of the strongly 
excited Yang-Mills field modes during the time when the 
field configuration is still far away from equilibrium and a 
quasi-classical description of the dynamics of the Yang- 
Mills field is appropriate, the ILEs are the most rele- 
vant Lyapunov exponents for the early thcrmalization at 
RHIC. 

Below we obtain the distribution of these three kinds of 
Lyapunov exponents. Since they govern the growth rate 
of the coarse grained entropy residing in the Yang-Mills 
field, they will allow us to estimate the equilibration time 
as r eq ~ AS/Sks, where A5* is the increase of entropy 
necessary for equilibration. 

Since classical CYM theory has no conformal anomaly 
(it does not know about Aqcd) all statistical quantities 
should scale like e"/ 4 , where e is the energy density and 
n is the mass dimension of that quantity. For example, 
the KS entropy has the mass dimension and scales as 
Sks oc e 1 / 4 . 

For the initial stage of high energy heavy ion collisions 
the relevant scale is the saturation scale Q s , which is 
related to the initial energy density in the color glass 
condensate (CGC) model as e = Qi/g 2 , implying that 
the time scale of very early dynamics is given by 1/Q S . 

Not surprisingly, Fries, Miiller and Schafer have indeed 
found that decoherence (which is probably the fastest 
mechanism for entropy production) happens indeed on 
this time scale @. 

However, they also found that decoherence can only 
generate a fraction of the entropy needed to justify a 
hydrodynamic treatment. 

The real-time gauge field dynamics discussed in this 
contribution is treated numerically introducing a spatial 
lattice with lattice constant a which accordingly has to 
be chosen as « ~ e -1 / 4 . We show that everything works 
out exactly in this manner. 

This paper is organized as follows. In Sec. |Hl we ex- 
plain the equations of motion in the CYM theory, and 
how we can obtain the eigenvalues of the Hessian in 
CYM. In Sec. IIII1 we show the eigenvalue distribution 
of the Hessian and its time evolution. Next we show the 
long-time evolution in terms of the maximum Lyapunov 
exponent. 



II. THEORETICAL BACKGROUND 

A. Chaotic dynamics of Yang-Mills fields 

In this section, following a brief review of previous re- 
sults, we discuss the method we use to analyze the com- 
plexity evolution in the classical Yang-Mills theory for 
an intermediate time duration. We first introduce the 
intermediate Lyapunov exponent which is applicable to 
general cases, and apply it to the classical Yang-Mills 
evolution. 

The chaotic properties of the classical evolution of 
Yang-Mills fields has been known and studied for a long 
time . Chaos was first observed in the infrared limit of 
the Yang- Mills theory Q ; later it was shown to exist also 
in the compact lattice version of the classical Yang-Mills 
theory [ltjj . The maximal global Lyapunov exponent may 
be related to the plasmon damping rate of the thermal 
pure Yang-Mills plasma [TTJ] . 

The global KS entropy of the compact lattice gauge 
theory (i.e. the rate of entropy growth close to thermal 
equilibrium) was shown to be extensive, i.e. proportional 
to the lattice volume [T2|, and the ergodic properties of 
the compact SU(2) lattice gauge theory were investigated 
numerically in detail by Bolte et al. [13J. 

Since we are here not interested in the ergodic prop- 
erties of the classical nonabelian gauge theory, but in its 
dynamical properties far off equilibrium, we will mostly 
make use of the non-compact formulation of the lattice 
gauge theory. In the following, we set the stage for our 
investigation by discussing three different kinds of insta- 
bility exponents, which capture different aspects of the 
dynamics of a nonlinear system with many degrees of 
freedom, such as the classical Yang-Mills field. 



B. Local and intermediate Lyapunov exponents 

For a simple "roll-over" transition, H = p 2 /2 — \ 2 x 2 /2, 
we have one positive and one negative Lyapunov ex- 
ponents, A and —A, which characterize both, the ki- 
netic instability and the entropy production. This is 
understood in the matrix form as follows. For a clas- 
sical trajectory, X(t) = (x(t),p(t)) T , we consider a sec- 
ond trajectory which differs a little in the initial con- 
dition. The equations of motion for the tangent vector 
SX(t) — (Sx(t),Sp(t)) T are written as, 

*M " (-1 J) g) • <«> 

where we have introduced short-hand notations, H x = 
dH/dx, H xp = d 2 H/dxdp, and so on. For an inverted 
harmonic oscillator, we put H xx = —A", H pp = 1, and 
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find, 



(6) 



This leads to an exponential expansion in the direction 
of Xx + p and an exponential contraction in the direction 
of —Ax + p. The entropy production rate in this simple 
case was analyzed by KMOS, who found to be given by 
dS/dt —> A for t co. 

In the case of many degrees of freedom, a similar struc- 
ture will appear as 

6X(t) = ( H £ x B„ \ sx{t) = H{t)sx{t) (7) 

Now the second derivatives should be regarded as ma- 
trices, e.g. (H xx )ij = d 2 H I 'dxidxj. We will refer to 
the matrix of second derivatives, Ti. as the Hessian in 
this paper. The eigenvalues A LLE of H are referred to as 
the local Lyapunov exponents (LLE). The LLE plays the 
role of a "temporally local" Lyapunov exponent, which 
specifies the departure of two trajectories in a short time 
period. 

If % is constant, i.e. in the absence of mode coupling, 
the LLEs arc identical with the Lyapunov exponents, and 
the KS entropy is defined as the sum of positive LLEs. In 
general, however, for a system with many degrees of free- 
dom, stable and unstable modes couple with each other. 
Thus, the LLE does not generally agree with the Lya- 
punov exponent in a long time period. In order to discuss 
the exponentially growing behavior of the fluctuation, we 
introduce the intermediate Lyapunov exponent (ILE). 

We can formally solve the equation of motion (JT)) for a 
finite time period At as, 



SX{t + At) = U(t, t + At)SX(t) 

pt+At 

U(t,t + At) =T cxp ' 



H(t + t')dt' 



(8) 
(9) 



where T denotes the time ordered product. Numerically, 
we can obtain the time-evolution operator U by the Trot- 
ter formula, 



U{t,t + At)=T J] U{t k -i,t k ) 

fe=l,JV 

-T \ ~ [l + n(t k -i)5i\ , (10) 



k=l.N 



where St = At/N. We diagonalize the time evolution 
matrix U and define the ILEs as, 

U D (t,t + At) =diag(e Al i LEAt ,e A " EAt ,...). (11) 

Liouville's theorem dictates that the determinant of the 
time evolution matrix U is unity, and thus the sum of all 
positive and negative ILEs is zero. After a long enough 
time for thcrmalization, the distribution of the ILEs is 



expected to converge to that of the global Lyapunov ex- 
ponents (GLE): 



ILE JA LLE (At^O), 
1 A GLE (At -)• oo) , 



(12) 



In general all three types of Lyapunov exponents, LLE, 
ILE, and GLE, yield different results. Here we are inter- 
ested in the rapid growth of the coarse grained entropy 
when the gauge field configuration is still far from equi- 
librium, but has already had sufficient time to sample a 
significant fraction of phase space. Our goal is not to cal- 
culate how the entropy grows when a configuration close 
to equilibrium relaxes further; this can be calculated re- 
liably in thermal quantum field theory. Instead, we focus 
below on the ILEs, and estimate the KS entropy as 



dS 



-£-S K8 - ^ A, 



ILE 



(13) 
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C. Classical Yang-Mills equation 

We consider the pure Yang-Mills theory in the tem- 
poral gauge, which permits a Hamiltonian formulation. 
The continuum Hamiltonian is given in terms of the phys- 
ical chromoelectric and chromomagnetic fields, Ef and 

B*=e ijk F? k ,by 



H 



J J \a,i a,i,j J 



(14) 



We now define the dimensionless variables on the lattice 
with lattice spacing a as (omitting vector and color in- 
dices) 



(15) 



A L = aA, E L = a 2 E, F L = a 2 F. 

The time variable is rescaled as 

t L = t/a. (16) 

The lattice spacing a is thus scaled out, and the dimen- 
sionless lattice Hamiltonian is defined as 



H L 



ag 2 H. 



(17) 



Here we make use of the fact that a rescaling of the 
Hamiltonian (by g 2 ) does not affect the classical equa- 
tions of motion. In the following we omit the superscript 
"L" . The Hamiltonian on the lattice is 



h =\T. e ^ x ) 2 + \ E w 2 . ( 18 ) 

x,a,i x,a,i,j 

F*(x) = diA^x) ~ djAfix) + £ r bc At(x)A%x) , 

(19) 



b.r 
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where di is the central difference operator in the i- 
direction, i.e., diA(x) = {A(x + i) — A(x — i)}/2. 
The classical equations of motion are given as, 

A?(x) = E?(x) , (20) 
Ef(x) = ]T OiW + £ f abc A)(x)F^(x) . (21) 

3 b,c,j 

There are two conserved quantities; total energy and 
charge. In the numerical simulation, we check that the 
total energy is strictly conserved along the time evo- 
lution. The charge conservation is expressed by non- 
Abelian Gauss' law, 

]T 8iE?(x) + ]T f abc A*(x)E?(x) = (22) 



where e L = (i? L ) / L 3 is the energy per site, i.e., the 
energy density in lattice units. 

A classical Yang-Mills theory on the lattice is a classi- 
cal system of 2L 3 (N 2 — 1) oscillators and has the thermal 
energy density 



= 2{N 2 C - 1)C L T L , 



(31) 



where Cl = ^2 k / L 3 is a numerical coefficient of order 
unity. The physical energy density of the lattice theory 
is 



£d(T) 



2(7V C 2 - 1)O l 



T 



(32) 



In the non-compact formalism, the lattice discritization 
violates the charge conservation in the magnitude of the 
field amplitudes. 

The Hessian of CYM theory is written as 



nj _ ( HeA HeE \ 

\-H AA -H AE J 
where the matrix elements arc 



H ee = S ab fiij8 x ,y , 
Hea = Hae = , 



(23) 



(24) 
(25) 



H A a = \5 ab P f ahC ® C + J2 f acd f bce R de i (26) 

c cde 

with 

P ~ V*x+l,y+3 ^x+i,y-~j ^x-i,y+j ^x-i,y-j) 

fe 

Q c =A1{y){5 XjV+ - - 8 XtU _ 3 ) - A%x){5 x ^ v - S x _ ly ) 

+ 5 l3 Y,{A%{x) + At{y)}{5 x+Ly - 5 x „ Ly ) 
k 

+ 2Ff j (x)6 x , y (28) 
R de ={-Al(x)A<*(x) + <% ]T A d k (x)At(x)}S x , y . (29) 

k 

On the L 3 lattice, the number of the eigenvalues is 6(N 2 — 
1)L 3 . 

D. Physical scale 

In order to fix the scale of the theory, we consider a 
physical volume V = a 3 L 3 in which the gauge field is 
thermalized at temperature T. The total energy is given 
by: 



{H) = Ve(T) 



(H L ) _ I?_ 

2 2 

9 a g z a 



(30) 



where we have used the relation T L = ag 2 T. On the 
other hand, the energy density in the weakly interacting 
thermal quantum Yang-Mills theory is 



e{T) = 2(N 2 - 1) 



d 3 k 



(2tt) 3 el k l/ T - 1 



2(7V C 2 -1)-T 4 



(33) 



The classical theory only applies to those modes of the 
continuum theory which are highly occupied and for 
which the quantum corrections are not too large. This 
condition imposes a lower limit on the lattice spacing of 
the classical theory. One can either argue that the two 
expressions for the energy density should coincide, or that 
a is the screening length of the corresponding quantum 
field theory. In both cases this leads to the relation 



a > 



T 



(34) 



where 9 is a numerical constant of order unity and T ~ 
e 1 / 4 is introduced as a measure of the energy density. 
(For example, e d = e leads to a = (30C L /7r 2 ) 1 / 3 /T ~ 
1.45/T.) 

The KS entropy growth rate, i.e. the sum of all positive 
Lyapunov exponents, in lattice units is given by 

4 L S = cks£ 3 (£ L ) 1/4 • (35) 
The KS entropy density in lattice units is thus 

^s=CKs(e L ) 1/4 . (36) 

The equilibrium entropy density of the classical Yang- 
Mills theory in the continuum with ultraviolet cut-off, is 
according to 



4e 



Se q (T) = -^ = 2(N^-l) 



4 Or. 



the same result in lattice units is 



fl W = a 3 Soq(T ) = 2(N 2 - 1) 



4C L 



(37) 



(38) 
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FIG. 1: Time evolution of the distance in SU(2) simulation FIG. 2: Time evolution in SU(2) simulation on 4 3 ,10 3 , and 
on 4 3 lattice. All scales are given in the lattice unit. 20 3 lattices with the same energy density. 



The equilibration time in lattice units is thus 

ri q L > = 4| = 2(N? 1)1^( £ L )-V4 (39) 
41 3CKS 

Finally, the physical equilibration time is 



= 2(iV c 2 - IjJgj^-V* . (40) 

III. CLASSICAL YANG-MILLS EVOLUTION 

A. Lyapunov exponents 

We first discuss the Lyapunov exponents obtained by 
the numerical simulations of SU(2) CYM systems. Ini- 
tial conditions are prepared with Ef (x) = and random 
Af(x) ^ 0. To see the chaotic time evolution, we mea- 
sured the "distances" between two gauge configurations: 



y x a,i a,i 

d ff = /EiE^^-E^w 2 } 2 ■ ( 42 ) 

The two gauge configurations are set to be very close to 
each other at the initial time t = 0. 

In Fig. [1] we show the numerical results on a 4 3 lat- 
tice. The energy density is e = 0.014. After a short 
time, the distance of two trajectories start to deviate, 
and exponentially grows in the intermediate time region 
(50 < t < 120). Later it saturates to a maximum value 
(t > 120). The exponential growth rate of the distance, 
i.e., the linear slope of hi Dff, in the intermediate time 
region is A^> ~ 0.04. This growth rate Ap is governed by 



the maximum Lyapunov exponent for a finite time pe- 
riod. In Fig. [2J we show the lattice size dependence of 
the time evolution. Apart from the irrelevant constant, 
the time evolution is almost insensitive to the lattice size. 
This is consistent with the expectation that the present 
lattice calculation simulates a piece of hot matter occu- 
pying a much larger volume. 

We calculated the ILEs by using the Trotter formula 
(fT0| . In the practical calculation, we have adopted the 
following expression, 

which contains an 0(St 2 ) term and coincides with l+WSt 
up to 0(5t). The determinant of this matrix is equal to 
unity and thus protects the symplectic property of the 
evolution. The eigenvalues are real or pure imaginary. 
These eigenmodes correspond to the exponentially grow- 
ing or damping mode and the oscillating mode, respec- 
tively. Since Liouville's theorem ensures the sum of the 
ILEs is zero, the positive and negative ILEs should ap- 
pear in a pairwise manner. 

We show the ILE distribution in Fig. [3] The gauge 
configuration is the same as in Fig. [Q The distribution 
at t = corresponds to the LLEs of the initial condi- 
tion. A positive (negative) LLE corresponds to the tem- 
porally local negative (positive) potential curvature, and 
the maximum LLE is larger than Aj> Within a short 
time period (0 < t < 5), the maximal ILE rapidly de- 
creases and the number of positive ILEs increases. As 
the distribution of ILEs no longer evolves for t > 50, 
the KS entropy is, therefore, also constant for t > 50. 
In this time region, the maximum ILE is A^f ~ 0.04, 
which is close to Aj> This fact means that the ILE does 
correspond to the growth rate for a finite time period. 
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FIG. 3: Distribution of the intermediate Lyapunov exponents. The total number of the eigenvalues is 1152. The right panel 
is a closeup of the largest 300 eigenvalues. 
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FIG. 4: The spectrum of gauge fields for different 

times. 



B. Growth of low- and high-momentum amplitudes 

Since the classical lattice theory is not ultra-violet 
(UV) safe, the energy is exhausted in this limit mostly 
by UV modes, which are sensitive to the lattice cutoff. 
We note that the classical theory at nonzero temperature 
has no well-defined continuum limit; e.g., the Raylcigh- 
Jeans formula gives an energy density that diverges in 
that limit. To wit, the thermal classical Yang-Mills the- 
ory on a lattice has no well defined continuum limit and 
the choice of lattice constant has physical significance, as 
discussed above. 

In order to examine further the role of UV modes we 
discuss next the momentum spectra of the gauge field 
A{p). We performed SU(2) simulations on a Nf = 16 3 
lattice with the energy density e = 0.014, which is the 
same setup as that in Sec. MI Al The distance Dff then 
exhibits a similar behavior as Fig. [TJ 

The gauge-field's spectra A(p) are obtained with 3-dim 
Fourier transformation of A(x), where momenta Ppijx = 
1,2,3) range from -(iV s /2-l) to N s /2 times 2ir/(N s a). 



(Note that due to our definition of the Laplacian, which 
extends to x ± 2a we get a factor 2ir/N s rather than 
tt/JV.). 

We average over direction and color, and 

show the time evolution of the spectrum A(p) of A(x) in 
Fig. |4l where the spectra are plotted as functions 

of v sin 2 pi + sin 2 p2 + sin 2 p%. Due to discretization of 
space one encounters doublers, which is why only half of 
the Brillouin zone is plotted in Fig. 

At t = 0, A(x) is randomly distributed and hence the 
A(p) are almost independent of \p\. After a short time 
the IR modes are strongly excited, and they dominate 
the exponential growth of the distance between gauge 
field configurations, as was expected. 

At low momenta our results approach the classical 
equilibrium (equipartition) distribution explains the ten- 
dency of our results at low momenta, but it is not com- 
pletely reached even in the IR modes at the last stage of 
the exponential growth (t ~ 150). 

Our results show, in addition, that at earlier times 
(t < 50) one has IR modes with very rapid growth, which 
appear to be the modes associated with the largest Lya- 
punov exponents. This would fit the usual assumption of 
a bottom- up thermalization [l9j , except that it is rather a 
pre-thermalization, because phase space is filled rapidly, 
but the full approach towards equilibrium sets in only 
with the linear phase, i.e. for t > 50. 

C. Equilibration time of SU(3) Yang-Mills theory 

Next, we discuss the SU(3) CYM theory. In Fig. \b\ 
we show the time evolution of Dff in SU(3) simulation 
on a 4 3 lattice for several energy densities. By changing 
the initial amplitude of Af(x), we calculated time evolu- 
tions with different energy densities. In Fig. [6j we show 
the ILE distributions after a long time period, which no 
longer change along time evolution. Only the positive 
eigenvalues are shown. These are qualitatively the same 
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FIG. 5: Time evolution in SU(3) simulation on a 4 3 lattice. 



as the SU(2) simulations. 

In Table Hand Fig. [3 we show the SU(3) results of the 
Lyapunov exponents: the exponential growth of the dis- 
tance A_d, the largest LLE Aj^, the sum of the positive 
LLEs A^", the largest ILE Aj^, and the sum of the 
positive ILEs Ag„ E . As discussed in the previous section, 
the Lyapunov exponents should scale as e 1 / 4 because of 
the conformal invariance. As shown in Fig. [3 Aj^ and 
A^f are indeed proportional to e 1 / 4 . Other Lyapunov 
exponents slightly deviate from this scaling. This is be- 
cause the change of the field amplitude is not exactly the 
conformal transformation, e.g., the dimcnsionless ratio of 
the electric energy density to the magnetic energy density 
is changed. Numerically, however, the best-fit prcf actor 
is not much affected by the difference of the exponent in 
the following accuracy. 

After all, we extract the Lyapunov exponent as a func- 
tion of the energy density from this approximate scaling. 
By fitting, we find that the numerical prefactors are 



\d 


-O.lXE 1 / 4 , 


(44) 


\LLE 

^max 


~ 1 x e 1 / 4 , 


(45) 


\LLE / r 3 
A sum / ^ 


~ 3 x e 1 / 4 , 


(46) 


A ILE 

max 


~ 0.2 x e 1 / 4 , 


(47) 


\ILE / r 3 
A sum/ 


~ 2 x e 1 / 4 . 


(48) 



Thus, the KS entropy density is 

_ xILE / r 3 ^ o .1/4 _ „ v 1/4 
SKS — A sum /-L — 1 x £ ~ C KS X £ 



(49) 



From this result, we can evaluate the equilibration time 
of the SU(3) CYM theory. We take 6 ~ {30C L /n 2 ) 1/3 ^ 
1.45 and C* L ~ 1 typical case. Then we get T = 
ag 2 T = ^(SO/tt 2 ) 1 / 3 ~ 6 and e L = 2(Y 2 -1)C L T L ~ 90. 
Inserting these numbers into Eq. (|40l) . we obtain 



5 
T 



^dclay 



(50) 



Here Td e i ay was introduced to take the initial phase into 
account, in which Dpp is more or less constant, because 



FIG. 6: Distribution of the intermediate Lyapunov exponents 
in SU(3) simulations. 
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FIG. 7: The SU(3) results of the Lyapunov exponents, Ad, 
A^x, A™, A££, and A££. The broken line is e 1/4 . 



the strongly growing modes are not yet relevant. One 
would expect that r^iay fulfills approximately [l6j 



6(Y 2 - 1)L 3 



3 Amax Tdeiay 



(51) 



which is indeed in good agreement with our numerical 
findings. When T ~ 350 McV, the equilibration time is 
r oq ~ 3 fm/c, with a systematic uncertainty which can 
easily account for a factor of two. 



TABLE I: The Lyapunov exponents in the SU(3) CYM the- 
ory. 



L 3 


e 


A D 


\ LLE 

^max 


\ LLE 
^Eum 


\ILE 
"max 


\ ILE 
''sum 


4 3 


0.05 


0.05 


0.55 


80 


0.06 


32 


4 :i 


0.20 


0.08 


0.77 


114 


0.11 


62 


4 :! 


0.86 


0.14 


1.14 


174 


0.20 


115 


4 :i 


3.16 


0.23 


1.64 


265 


0.32 


191 


4 3 


19.4 


0.39 


2.89 


474 


0.59 


328 


4 :i 


90.9 


0.59 


4.60 


708 


0.95 


474 



8 



If entropy is produced by very rapid processes, especially 
by decoherence, before the non-linear dynamics analysed 
by us reaches the phase of linear entropy growth, the 
real thermalization time is correspondingly shorter. In 
@ this decoherence entropy was estimated to be roughly 
1/3 of what is needed by the hydrodynamical initial con 



ditions. This is consistent with results obtained in [21 1 
in kj_ factorized perturbation theory. In that calculation 
the full observed particle number at central rapidities is 
only reached from decohering the Color Glass Conden- 
sate after introducing a normalization factor. Without 
that factor one would obtain between one half and one 
quarter of the total particle number. On the other hand, 
simulations of the combined decoherence and non-linear 
dynamics stage of the glasma in a longitudinally expand- 
ing, boost-invariant geometry reaches 80% of the final 
particle number [22j ■ All of this indicates that it is prob- 
ably a good guess to assume that non-linear gauge field 
dynamics has to generate about 2/3 of the entropy re- 
quired by thermal equilibrium and that the thermaliza- 
tion time is thus rather of the order of 2 fm/c. 



IV. SUMMARY 

The main aim of this paper is to understand the fast 
thermalization deduced for high energy heavy ion collis- 
sions, which is, in fact, one of the largest unexplained 
puzzles in RHIC physics. We argue that entropy genera- 
tion plays the key role in this context. 

The overall picture of entropy generation in heavy ion 
collins is involved. While some part of the entropy is 
produced from the decoherence at very early times, i.e. 
times of order 1/Q S , see [f|, most of entropy required by 
the initial condition for the hydrodynamic phase must be 
generated within the first fm/c by noncquilibrium gluon 
dynamics. Entropy generation of quantum systems al- 
ways requires coarse graining. Coarse-graining in turn 
relates it to the Lyapunov exponents, see As the lat- 
ter can be studied in the corresponding classical gauge 
theories so can entropy production in total. 

More precisely, the entropy production rate in classi- 
cal systems is given by the Kolmogorov-Sina'i (KS) en- 
tropy, defined as the sum of positive Lyapunov expo- 
nents. (The KS entropy describes the entropy produc- 
tion also in quantum systems when the coarse graining 
is introduced with a minimum wave packet Q.) 

We have investigated classical Yang-Mills (CYM) dy- 
namics in the noncompact (A, E) scheme. We started 
from random initial conditions and studied the result- 
ing spectrum of Lyapunov exponents. We found that 
their properties change with time in a characteristic man- 
ner and identified three distinct regimes: A short time 
regime, in which the system has not yet sampled a large 
fraction of phase space, a late time regime in which the 
system is already close to thermal equilibrium and has 
sampled basically all of phase space, and an intermedi- 
ate regime which is dominated by non-linear gauge field 



dynamics. 

We have developed a method, making use of Trot- 
ter formula, to evaluate the Lyapunov exponent in the 
intermediate time scale (intermediate Lyapunov expo- 
nent; ILE), which is the relevant time scale for the 
problem of thermalization in heavy ion collisions, and 
determined the entropy production rate (Kolmogorov- 
Sina'i entropy) . The obtained equilibration time scales as 
r cq cx e -1 / 4 oc 1 /T + Tdday , where £ is the energy density 
and Tdeiay is the typical time to reach the intermediate 
time range, which we also determined. In total the ther- 
malization time is around 2 fm/c for T = 350 MeV, if 
one assumes that 1/3 of the required entropy is gener- 
ated by decoherence, with rather substantial systematic 
uncertainties. The most important source of uncertainty 
is related to the choice of lattice constant a. Since CYM 
has conformal symmetry, the physical scale setting is pro- 
vided by the discretization scale a which thus acquires 
physical importance. The choice of a is not free of ambi- 
guities. Different arguments all lead to the form a = ce 1 / 4 
with a constant c of the order of one but varying within 
a factor of two. 

One also finds that the e dependence of a is crucial 
to obtain the correct power scaling for all quantities of 
interest from CYM. 

In the course of these investigations it was crucial to 
understand the qualitative differences between the dif- 
ferent time ranges and corresponding Lyapunov spectra, 
which also allows to reconcile previously not understood 
observations [20| . 

A thermalization time of roughly 2 fm/c is some- 
what larger than the phcnomenologically preferred value. 
However, the inclusion of quarks will reduce this number 
and could well bring it into the phcnomenologically pre- 
ferred range. In addition strong electric field in the initial 
condition together with the magnetic field, and longitu- 
dinal (Bjorken) expansion may promote faster equilibra- 
tion. 
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